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Preface 


This  study  began  with  the  investigation  of  oscilla¬ 
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Abstract 


The  two  dimensional  transient  heat  (diffusion) 
equation  with  Dirichlet  boundary  conditions  was  solved 
using  the  Duf ort-Frankel ,  Saul'ev,  and,  Exponential 
(Power-law)  finite  difference  schemes.  All  methods  were 
investigated  for  oscillatory  behavior  and  comparisons  of 
accuracy  made. 

To  predict  the  time  step  at  which  oscillatory 
behavior  would  occur,  the  coefficient,  matrix,  and  probabil¬ 
istic  methods  of  stability  analysis  were  utilized.  At 
time  steps  greater  than  the  square  of  the  mesh  divided  by 
the  thermal  diffusivity,  oscillatory  solutions  were  appar¬ 
ent  in  both  the  Duf ort-Frankel  and  Saul’ev  schemes.  The 
exponential  method,  as  predicted,  did  not  oscillate  for  any 
size  time  step. 

Although  the  exponential  scheme  was  the  most  accu¬ 
rate  at  large  time  steps,  the  solution  still  contained 
enough  error  to  be  unusable  in  many  engineering  applica¬ 
tions.  At  small  time  steps,  all  methods  were  more  accu¬ 
rate  than  the  fully  implicit  formulation. 

The  exponential  method  was  found  to  be  the  slowest 
computationally.  The  Saul'ev  scheme  proved  to  be  the 
fastest  while  still  achieving  the  required  degree  of 
accuracy . 
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OSCILLATORY  BEHAVIOR  OF  FINITE  DIFFERENCE  METHODS 
FOR  THE  SOLUTION  OF  THE  TWO  DIMENSIONAL 
TRANSIENT  HEAT  (DIFFUSION)  EQUATION 


I.  Introduction 


Background 

Because  only  the  simplest  of  boundary  value  prob¬ 
lems  can  be  solved  by  exact  analytical  methods,  the  engi¬ 
neer  must  rely  on  finite  difference  methods  in  solving 
"real  world"  problems  with  complicated  and  irregular 
geometries.  Finite  difference  approximations  allow  a  dif¬ 
ferential  equation,  together  with  its  associated  boundary 
and  initial  conditions,  to  be  transformed  into  a  set  of 
algebraic  equations  readily  solved  utilizing  a  high  speed 
digital  computer. 

Clearly,  for  finite  difference  approximations  to 
be  of  practical  use,  a  high  degree  of  accuracy  must  be 
demanded.  It  is  well  known  that  difference  approximations 
for  derivatives  contain  "truncation"  errors  arising  from 
the  neglect  of  the  higher  order  terms  in  its  Taylor  series 
representation.  "Round-off"  errors  are  also  introduced 
since  only  a  finite  number  of  digits  can  be  stored  in  com¬ 


puter  memory. 


A  finite  difference  scheme  is  said  to  be  "stable" 
if,  as  time  progresses,  these  errors  decay  rather  than 
grow  in  an  unbounded  fashion  (2:98).  "The  inexperienced 
user  is  often  surprised  when  stable  methods  produce 
physically  unrealistic  solutions  that  oscillate  with  time" 
(9:27).  Although  in  stable  methods,  these  oscillations 
decrease  with  time,  they  may  produce  sufficient  inaccura¬ 
cies  to  destroy  the  solution  (9:27).  Therefore,  a  thorough 
understanding  as  well  as  a  means  of  predicting  oscillatory 
behavior  is  essential  in  employing  a  finite  difference 
method  with  any  degree  of  confidence. 

The  Two  Dimensional  Transient 
Heat  (Diffusion)  Equation 

Of  particular  interest  to  the  nuclear  engineer  are 
boundary  value  problems  involving  the  transient  heat 
(diffusion)  equation.  The  equation  in  two  dimensional 
form  is: 


where  a  is  thermal  diffusivity  and  T  is  temperature. 
Problem  Statement 

The  objective  of  this  thesis  is  to  determine  if 
numerical  oscillations  exist  in  the  following  finite 


difference  approximations  to  the  two  dimensional  transient 
heat  (diffusion)  equation: 

1.  Dufort-Frankel 

2 .  Exponential 

3 .  Saul ' ev 

In  addition  to  the  above  methods,  solutions  found 
by  the  pure  implicit  method  are  computed  for  purposes  of 
accuracy  comparison. 

Scope 

The  boundary  value  problems  investigated  are 
limited  to  the  two  dimensional  transient  heat  equation 
applied  on  a  unit  square.  Dirichlet  boundary  conditions 
and  a  thermal  diffusivity  of  one  are  assumed.  The  unit 
square  is  partitioned  into  a  grid  of  81  interior  nodes 
with  ax  »  Ay-  The  analysis  of  oscillatory  behavior  is 
conducted  by  following  the  temperatures  of  two  interior 
nodes.  All  nodes  will  be  included  in  the  error  analysis. 

Approach 

A  previous  thesis  by  Kropf  demonstrated  that  stable 
oscillations  existed  in  the  Alternating  Direction  Implicit 
Method  of  Peaceman  and  Rachford  when  applied  to  a  unit 
square  (6:56).  After  confirming  these  results,  a  litera¬ 
ture  search  was  conducted  to  find  other  difference  methods 
that  could  be  subject  to  oscillatory  behavior.  The  Dufort- 
Frankel,  exponential,  and  Saul'ev  finite  difference  schemes 
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were  selected  and  tested  for  oscillations  using  the  coeffi¬ 
cient,  probabilistic,  and  matrix  methods.  Two  boundary 
value  problems  were  then  numerically  solved  using  the  above 
finite  difference  methods  and  investigated  for  the  oscilla¬ 
tions  predicted. 

Sequence  of  Presentation 

General  algebraic  formulas  for  the  finite  differ¬ 
ence  methods  considered  are  presented  in  Chapter  II.  The 
coefficient,  probabilistic,  and  matrix  methods  of  predict¬ 
ing  oscillations  are  introduced  in  Chapter  III. 

Chapter  IV  presents  numerical  solutions  for  the 
two  problems  as  well  as  accuracy  comparisons  for  all 
methods.  Finally,  conclusions  and  recommendations  are 
offered  in  Chapter  V. 


II.  Difference  Approximations  to  the  Transient 
Heat  (Diffusion)  Equation 

The  Dufort-Frankel,  Saul'ev,  and  exponential  finite 
difference  approximations  will  be  investigated  for  oscilla¬ 
tory  behavior  and  solution  accuracy.  The  Dufort-Frankel 
and  Saul'ev  difference  schemes  are  explicit  methods  that 
are  known  to  be  unconditionally  stable  in  two  dimensions 
(7:237).  A  relatively  new  implicit  method,  the  exponential 
scheme  of  Patankar  and  Baliga,  is  also  considered.  In  the 
one  dimensional  case,  this  method  has  been  shown  to  be  free 
of  numerical  oscillations  and  possess  superior  accuracy  for 
any  size  time  step  (9:35).  Here,  the  analysis  is  extended 
to  two  dimensions. 

Standard  difference  notation  is  used  in  presenting 
alj.  methods.  Namely,  T(x,y,t)  =  T (iAx,  jAy ,nAt)  =  T(i,j,n). 
Further,  the  following  nomenclature  is  adopted: 

H  =  AX  =  Ay 

and  R  =  aAt/H2  =  At/H2 

Dufort-Frankel  Explicit 
Approximation 

The  Dufort-Frankel  approximation  was  developed  to 
circumvent  the  inherent  instability  of  Richardson's  three 
level  explicit  method  (7:238).  The  formula,  for  two  dimen¬ 
sions,  is  (App  A) : 
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T(i, j,n+l) (1+4R)  =  2R  [T (i+1 , j ,n)  +  T(i-l,j,n) 
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+  T(i, j+l,n)  +  T  ( i, j-1 , n) ]  (1) 

+  (1-4R)  T (i, j , n-1) 

The  method  is  effectively  explicit  if  temperatures 
at  both  the  (n-1)  and  (n)  time  steps  are  known  or  in  the 
computationally  limiting  case  of  R  equal  to  1/4  (7:239). 
Given,  as  an  initial  condition,  the  temperatures  at  the 
(n-1)  time  level,  an  alternate  method  is  needed  to  advance 
to  time  level  (n)  (7:239). 

Saul 1 ev  Explicit  Approx ima t ion 

The  Saul'ev  scheme,  like  the  Dufort-Frankel,  is 
explicit  under  certain  conditions.  However,  the  Saul'ev 
scheme  employs  alternate  formulas  for  each  successive  time 
step.  To  advance  to  the  (n+1)  time  level  the  equations  are 
(7:237) : 

T (i, j , n+1)  (1+2R)  =  R  (  T ( i, j+1 ,n)  +  T(i+l,j,n) 

+  T(i-1, j ,n+l)  +  T (i, j-1, n+1) ]  (2a) 

+  (1-2R)  T(i, j ,n) 

To  proceed  from  the  (n+1)  to  the  (n+2)  time  step, 
the  equations  are  (7:238): 

T(i,j,n+2)  (1+2R)  =  R  (  T(i+1, j ,n+2)  +  T(i,j+l,n+2) 

+  T ( i-1 , j , n+1)  +  T(i, j-1, n+1) ]  (2b) 

+  (1-2R)  T ( i , j , n+ 1 ) 
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Equation  (2a)  becomes  explicit  when  T(i-l,j,n+l) 
and  T(i,j-l,n+l)  are  known.  If  using  rectangular  geometry, 
this  is  accomplished  by  initiating  the  calculations  in  the 
lower  left  corner  and  proceeding  up  the  grid  from  left  to 
right  (7:195).  Similarly,  Equation  (2b)  can  be  "made1' 
explicit  by  starting  in  the  upper  right  corner  and  proceed¬ 
ing  down  the  grid  from  right  to  left  (7:195). 

Exponential  Approximation 

The  exponential  method  of  Patankar  and  Baliga  was 
shown  to  be  non-oscillatory  in  the  one  dimensional  case 
(9:37).  A  thorough  development  of  the  method  can  be  found 
in  reference  9.  The  difference  formula  in  two  dimensions  is 

T(i,j,n+1)  ( l+4Rf )  -  Rf [T(i+1, j ,n+l)  +  T(i-l,j,n+l) 

+  T(i, j+l,n+l)  +  T(i, j-l,n+l) ] 

=  R(l-f) [T(i+1, j,n)  +  T ( i-1 , j , n)  (3) 

+  T{i, j+1 , n)  +  T ( i, j-1 , n) ] 

♦  [1-4R  (1-f ) ]  T(i, j ,n) 

where 

f  =  (1/  [  1-exp (-4R) ] )  -  (1/ (4R) ) 

To  achieve  a  computationally  faster  algorithm,  the 
authors  replaced  the  exponential  term  with  the  following 
power  law  approximation  (9:32): 


For  0  <  4R  <  10  f  *  1  -  [l/(4R)j  [ 1- ( 1-0 . 1 (4R) ) 5 J 

For  4R  >  10  f  »  1  -  [ 1/ ( 4 R) ] 
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This  difference  scheme,  unlike  the  others,  results 
in  a  system  of  equations  to  be  solved  simultaneously. 

Methods  of  Solution 

The  Dufort-Frankel  and  Saul'ev  explicit  methods  are 
solved  in  a  straightforward  manner  with  one  equation  for 
each  unknown.  However,  both  methods  are  explicit  under 
certain  conditions. 

To  "bootstrap"  the  three-level  Dufort-Frankel 
scheme  through  the  first  two  time  steps,  another  finite 
difference  method  is  needed.  The  "pure  implicit"  method 
was  chosen  because  of  its  characteristic  non-oscillatory 
behavior.  With  the  first  two  time  steps  known,  the  method 
subsequently  became  explicit. 

For  an  explicit  Saul’ev  scheme  to  be  realized,  the 
order  in  which  the  nodal  points  are  calculated  is  impor¬ 
tant.  For  each  odd  time  step,  equation  (2a)  is  used, 
starting  in  the  lower  left  corner  and  proceeding  up  the 
grid  from  left  to  right.  Alternatively,  Equation  (2b)  is 
used  for  each  even  time  step  starting  in  the  upper  right 
corner  and  moving  down  the  grid  from  right  to  left. 

The  exponential  method  resulted  in  a  series  of 
equations  to  be  solved  simultaneously.  These  were  of  the 
form 
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A  x  =  b 
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where  x  is  a  column  matrix  of  temperatures  at  time  (n+1) 
and  b  is  a  column  matrix  of  the  associated  boundary  condi¬ 
tions  and  temperatures  at  time  (n) .  For  N  interior  nodes, 
a  real,  symmetric,  N  x  N  matrix  results  for  A.  The  most 
computationally  efficient  method  of  solving  this  matrix 
equation  is  direct  Gaussian  elimination  (13:105). 
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III.  Oscillatory  Behavior  and  Accuracy  of  Methods 

Unfortunately,  the  stability  of  a  finite  difference 
method  does  not  guarantee  solution  accuracy  or  even  physical 
realism  (9:27).  Stability  merely  insures  that  numerically 
induced  errors  will  not  grow  in  successive  levels  of  compu¬ 
tation  (2:98).  For  computational  efficiency,  the  largest 
time  step,  which  produces  the  required  level  of  accuracy, 
should  be  utilized.  If  too  large  a  time  increment  is  used, 
numerical  oscillations  may  result,  destroying  solution 
accuracy  (9:27).  Attention  is  now  turned  to  three  methods 
of  stability  analysis  that  are  of  use  in  predicting  oscilla¬ 
tory  behavior. 


o- 

«a 


Coefficient  Method 

The  method  presented  here  was  proposed  in  a  pre¬ 
vious  paper  (6:10).  It  is  a  direct  extension  of  the  one 
dimensional  case  offered  by  Myers  (8:282).  By  applying 
Equations  (1),  (2),  and  (3)  to  a  one  nodal  point  grid  with 
homogeneous  boundary  conditions,  the  following  ratios  are 
obtained: 

Dufort-Frankel 

T(i,j,n+1)/T(i,j,n-1)  =  ( 1-4R) / ( 1+4R) 


i 

S 
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Saul ' ev 

T(i,j,n+1)/T(i,j,n)  =  (1-2R) / (1+2R) 

10 


Exponential 


T(i,j,n+1)/T(i,j,n)  =  exp(-4R) 

The  value  of  this  ratio  (p)  will  reveal  the 
behavior  of  the  corresponding  finite  difference  scheme 
(8:282)  : 

For  0  <  p  <  1  No  oscillations,  stable  solution 

For  -1  <  p  <0  Stable  oscillations,  stable  solution 

For  p  <-l  Unstable  oscillations,  unstable  solution 

It  is  readily  apparent  that  all  ratios  lie  between 
-1  and  1,  and  therefore,  all  methods  are  unconditionally 
stable.  However,  a  region  of  stable  oscillations  is  pre¬ 
dicted  for  both  the  Duf ort-Frankel  and  Saul'ev  schemes 
(see  Figure  1) . 

For  the  Duf ort-Frankel  method,  oscillatory  behavior 
is  predicted  when 

R  >  1/4 

This  value  of  R  causes  the  ratio  to  become  negative  and 
stable  oscillations  should  be  observed. 

Similarly,  the  Saul'ev  scheme  is  predicted  to 


oscillate  when 


Note  the  method  predicts  only  smooth,  stable 
behavior  for  the  exponential  scheme  since  the  ratio  is 
always  positive  and  less  than  unity. 

It  is  important  to  remember  that  since  no  more 
than  one  node  and  homogeneous  boundary  conditions  were 
assumed,  the  coefficient  method  can  only  approximate  the 
values  of  R  that  cause  numerical  oscillations  (8:289). 

Matrix  Method 

Unlike  the  coefficient  method,  the  matrix  approach 
deals  with  the  more  general  case  in  which  there  is  more  than 
one  node.  Myers  presents  this  approach  employing  a  two 
node  grid  (8:286).  The  method  here  will  be  extended  to 
four  nodes. 

Consider  a  system  of  N  finite  difference  equations 
that  can  be  expressed  as 

u,n+1)  =  A  i,nl 

where  u^n+^  and  u^  are  column  matrices  comprised  of 
temperatures  at  times  (n+1)  and  (n)  respectively.  The 
spectral  radius  of  A  then  determines  the  stability  of  the 
system  (11:84) .  Namely,  if  the  eigenvalue  of  largest 
modulus  is  less  than  or  equal  to  unity,  the  system  will 
be  stable. 

For  example,  when  the  exponential  scheme  is  written 
in  the  above  form  and  a  four  node  grid  is  assumed,  the 


asp 


eigenvalues  of  A  are  found  to  be  (App  B) : 

A 1  =  A  2  =  exp (-4R) 

A  3  =  1- (6R/ (l+6Rf) ) 

A  4  =  1- (2R/ (l+2Rf ) ) 

It  is  easily  verified  that  A 4  is  the  dominant 
eigenvalue  and  lies  between  0  and  1.  Using  the  two  node 
model  of  the  Euler  method,  Myers  finds  the  value  of  R  which 
would  make  the  dominant  eigenvalue  negative,  and  states 
that  this  is  the  "limit  for  stable  oscillations"  (8:286). 

In  applying  this  reasoning  to  the  exponential  case,  one 
concludes  that  for  any  value  of  R  the  dominant  eigenvalue 
is  always  positive  and  less  than  unity.  Thus,  the  matrix 
approach  predicts  only  steady,  stable  numerical  behavior  in 
agreement  with  the  coefficient  method. 

It  is  easy  to  see  that  this  method  becomes  pro¬ 
hibitively  laborious  when  large  systems  of  equations  (nodes) 
are  considered  (8:286). 

Probabilistic  Method 

A  probabilistic  method  of  stability  analysis  is 
suggested  by  Kaplan  (5:459).  If  the  coefficients  of  an 
explicit  finite  difference  scheme  are  given  probabilistic 
interpretations,  it  has  been  shown  that  the  stability 
limit  can  be  derived  (5:459).  This  concept  can  be  taken 
a  step  further  by  utilizing  the  method  to  determine 
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oscillatory  behavior  in  explicit  methods.  The  method  will 
be  applied  to  the  Duf ort-Frankel  scheme  since  it  has  not 
been  demonstrated  to  be  of  use  in  implicit  methods. 

By  giving  the  coefficients  found  in  the  Dufort- 
Frankel  scheme  a  probabilistic  interpretation.  Equation  (1) 
can  be  rewritten 

T ( i, j , n+1)  =  P"  T ( i+1 , j , n)  +  P*  T(i-l,j,n) 

+  P~  T ( i, j+1 , n)  +  P*  T ( i, j-1 , n) 

+  PQ  T(i, j ,n~l) 

where 

Px  =  Px  =  Py  =  Py+  =  <2R)/<1+4R) 

P  =  (1-4RJ  /  (1+4RJ 
o 


i 


Moreover,  Kaplan  points  out  that  these  probability 
coefficients  must  satisfy  two  conditions  (4:460): 

1.  No  probability  or  coefficient  can  be  negative 

2.  The  sum  of  the  probabilities  or  coefficients 
cannot  exceed  unity 

By  simply  adding  the  coefficients  found  in  Equa¬ 
tion  (1),  condition  two  is  found  to  be  satisfied.  However, 
condition  one  implies  that  Pq  must  be  non-negative  which 
occurs  only  when 


R  <  1/4 


The  method  then  predicts  that  oscillations  may 
occur  when  R  is  greater  than  1/4  and  is  in  agreement  with 
the  coefficient  method. 


Accuracy  of  Methods 

Since  finite  difference  methods  are  "built"  by 

neglecting  the  higher  order  terms  of  its  Taylor  series 

representation,  truncation  error  is  inevitable.  The 

Duf ort-Frankel  scheme  possesses  truncation  error  of  order 
2  2  2 

H  +  At  +  (At/H)  (7:157).  An  improvement  in  the  trunca¬ 
tion  error  of  the  Saul'ev  method  can  be  realized  by  using 

Equations  (2a)  and  (2b)  in  ?n  alternating  fashion  (7:199). 

2 

If  each  is  used  alone,  the  error  is  of  order  (H  +  At  ) 


(7:199).  When  used  together  in  alternating  directions, 

2  2  2 

the  error  is  reduced  to  order  H  +  At  +  (At/H)  (7:199) 
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The  truncation  error  for  the  exponential  scheme  is 
not  listed  in  the  literature.  However,  an  approximation 
can  be  found  by  noting  that  for  certain  values  of  f.  Equa¬ 
tion  (3)  is  equivalent  to  other  conventional  finite  differ¬ 
ence  schemes  with  known  truncation  error.  With  f  equal  to 
0.5,  Equation  (3)  corresponds  to  the  Crank-Nicolson  method 
(9:29).  An  f  of  1.0  similarly  leads  to  the  fully  implicit 
formulation  (9:29).  Further,  it  is  easily  verified  that 
for  large  values  of  R,  f  approaches  1.0  in  value.  Con¬ 
versely,  as  R  is  made  arbitrarily  small,  f  approaches  0.5. 
Thus,  for  small  values  of  R,  the  exponential  scheme  should 
have  truncation  error  approaching  the  Crank-Nicholson 
(H2  +  At2)  (10:189)  and  at  large  R,  that  of  the  fully 
implicit  (H2  +  At)  (10:189). 

In  addition  to  truncation  error,  oscillatory 
behavior  can  cause  considerable  inaccuracy.  Substantial 
errors  are  then  predicted  for  both  the  Duf ort-Frankel  and 
Saul'ev  methods  for  large  values  of  R. 

Because  the  computer  can  retain  only  a  finite 
number  of  digits,  round-off  errors  are  introduced.  As  the 
number  of  computations  are  increased,  the  round-off  error 
may  grow,  possibly  nullifying  any  accuracy  advantage  gained 
by  refining  time  step  or  mesh  size  (6:15). 


Error  Analysis 

To  determine  the  accuracy  associated  with  each 
finite  difference  method,  maximum  and  root  mean  square 
errors  were  computed.  Root  mean  square  error  was  calcu¬ 
lated  according  to  (12:1023-1024)  . 

J  I 

Erms  =  2  ^[(Ta(i,j,n)  -  T  (i,  j  ,n)  )  2/IJ]  ** 

j-1  i-1 


where 

Ta(i,j,n)  =  analytical  solution  at  point  (i,j,n) 

T(i,j,n)  =  finite  difference  solution  at  point 
(i, j ,n) 

I  =  number  of  nodes  in  x  direction 
J  *  number  of  nodes  in  y  direction 

The  maximum  error  at  any  nodal  point  was  also  cal¬ 
culated.  This  gives  an  indication  of  overall  grid  accu¬ 
racy  and  is  computed  as 

Emax  *  max  |  Ta(i,j,n)  -  T(i,j,n)  | 
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IV.  Numerical  Results 


Numerical  solutions  obtained  by  the  Duf ort-Frankel , 
Saul'ev,  and  exponential  schemes  are  presented.  Two  sample 
problems,  each  with  a  known  analytical  solution,  are  solved 
using  each  method.  In  addition,  solutions  obtained  by  the 
fully  implicit  method  are  provided  for  comparison.  The 
sample  problems  and  their  analytical  solutions  are  first 
presented.  Next,  oscillatory  behavior  is  investigated  for 
all  methods.  Finally,  root  mean  square  (Erms)  and  maximum 
(Emax)  error  are  given  to  quantify  solution  accuracy. 


Problem  #1 

The  first  problem  and  its  analytical  solution  were 
successfully  used  in  a  previous  thesis  to  investigate 
oscillatory  behavior  of  the  Peaceman-Rachf ord  alternating 
direction  implicit  method  (6:18-19).  The  problem  is: 


3T(x,y,  t)  _  3 2T (x , y , t)  32T(x,y, t) 

2  2 
3 t  3x*  3y^ 


with 


T(x,y,0)  =  0 


T (x, 0 , t)  =  T (x, 1 , t)  =  T(l,y,t)  =  T(0,y,t)  =  400 


vV Vi ; m  'V 


R 
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0  <  x  <  1  and  0  <  y  <  1 
The  analytical  solution  is  (6:19) 
T  (x,y ,  t)  =  400  (1-f (x,y ,  t) ) 


where 


f(x,y,t)  =  g(x,t)  h (y , t) 


La  2  2 

—  exp(-n  it  t)  sin 


(mrx) 


Z4  2  2 

—  exp  (-n  it  t)  sin  (niry) 


n  *  If  3  f  5  ••• 


Problem  #2 

The  second  problem  considered  is: 


3T(x,y , t)  _  32T (x,y , t)  32T (x,y,t) 

2  2 

at  ax  ay 


with 


where 


T (x,y , 0)  =  0 

T(x,0,t)  =  T(0,y,t)  =  T(l,y,t)=  0;  T(x,l,t)  =  400 

0  <  x  <  1  and  0  <  y  <  1 


m 

5 mom 

This  problem  is  identical  to  one  found  in  Myers 


solved  for  the  steady  state  case  (8:127).  The  solution  for 
the  transient  case  is  found  to  be: 

T  (x ,  y ,  t)  =  f  (x,y, t)  +  g(x,y) 


where 

QO  00 

^  2  2  2  2 
f  (x,y,  t)  =  2^  2^  sin  (mux)  sin(nuy)  exp(-(n  it  +  m  rr  )  t) 

m=l  n=l 

n  =  1,  2/  3 ,  ... 
m  =  1/  2 ,  3 1  ... 


and 


g(x,y)  =  -  [(sin(lnx)  sinh  (l^y) )  /  (1  sinh(lir))] 

1=1 


1=1,  3,  5, 


B 


mn 


-rr 


g(x,y)  sin(n^y)  sin(m^x)  dx  dy 


The  steady  state  solution,  g(x,y),  is  found  in 
Myers  (8:129).  The  transient  portion  is  found  by  separa¬ 
tion  of  variables  and  superposition  (1:130). 


Oscillatory  Solutions 

Solutions  for  each  method  were  computed  using  a 
square  grid  with  H  =  ax  =  Ay  =  0.1  (Figure  2) .  To  graphi¬ 
cally  display  oscillatory  behavior,  solutions  for  nodes  57 
and  41  are  shown  for  selected  values  of  R  and  corresponding 
time  step  AT.  Temperatures  are  listed  in  degrees.  The 
specific  unit,  whether  °K  or  °C  will  obviously  follow  from 
the  choice  of  thermal  diffusivity  (a),  here  equal  to  one. 


Problem  #1 

As  predicted  in  Chapter  III,  oscillatory  solutions 
were  found  for  both  the  Dufort-Frankel  (DF)  and  Saul'ev  (SV) 
difference  methods  when  particular  values  of  R  were  util¬ 
ized.  Figures  3  and  4  demonstrate  the  oscillations  pro¬ 
duced  by  the  DF  for  R  =  1.0.  Although  the  oscillations 
decrease  with  time,  they  are  still  large  enough  (>  100 
degrees)  to  render  the  solution  unusable.  For  values  of 
R  below  0.25,  only  steady,  smooth  behavior  was  observed. 
Figures  5  and  6  show  the  DF  scheme  to  be  free  of  oscilla¬ 
tions  for  R  =  0.2. 

The  Saul'ev  approximation  was  found  to  be  free  of 
oscillations  for  R  =  1.0  (Figures  7  and  8).  However, 
oscillatory  behavior  became  evident  when  R  was  increased 
to  5.0  (Figure  9).  This  was  to  be  expected  since  the 
"threshold"  for  oscillations  in  the  SV  (R  =  0.5)  is  higher 
than  in  the  DF  method  (R  =  0.25). 
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FIG  3.  NUMERICAL  SOLUTIONS  PROBLEM  *1  (NODE  41) 
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FIG  4.  NUMERICAL  SOLUTIONS  PROBLEM  *1  ( NOOE  57 


FIG  8.  NUMERICAL  SOLUTIONS  PROBLEM  *1  (NODE  411 
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Oscillatory  behavior  could  not  be  induced  in  the 
exponential  scheme  (EXP)  no  matter  how  large  the  time  step. 
Figures  10  and  11  show  the  solution  for  R  =  1.0.  Solu¬ 
tions  for  R  =  10.0  (AT  =  0.1)  are  displayed  in  Figures  12 
and  13.  Again,  the  solutions  are  free  of  numerical  fluctua¬ 
tions. 

Problem  #2 

In  general,  problem  #2  results  confirmed  the  previ¬ 
ous  findings.  Namely,  that  both  the  DF  and  SV  methods  pro¬ 
duce  oscillatory  solutions  for  values  of  R  in  the  ranges 
predicted.  Numerical  oscillations  are  again  seen  for  the 
DF  method  at  R  =  1.0.  A  comparison  of  Figures  2  and  4 
with  Figures  14  and  15  reveal  similar  oscillations  that 
decrease  with  time. 

Unlike  before,  the  SV  scheme  now  displays  small 
oscillations  for  R  =  1.0  (Figures  16  and  17).  When  the 
value  of  R  is  increased  to  5.0,  the  oscillations  grow  in 
magnitude  and  gross  error  results.  However,  stability  is 
maintained  since  the  oscillations  decay  with  time  (Figure 
18)  . 

Once  again,  the  exponential  method  produced  only 
smooth,  steady  solutions  for  all  values  of  R.  Figures  19 
and  20  show  the  solution  to  be  non-oscillatory  for  R  =  1.0. 
Note  also  that  the  solution  is  generally  more  accurate  than 
the  implicit  method. 


FIG  11.  NUMERICAL  SOLUTIONS  PROBLEM  '1  ( NOOE  57 
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FIG  14.  NUMERICAL  SOLUTIONS  PROBLEM  *2  ( NOOE  41) 
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FIG  19.  NUMERICAL  SOLUTIONS  PROBLEM  *2  f  NO OE  57) 
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Accuracy  of  Solutions 

To  quantitatively  examine  the  effects  of  oscilla¬ 
tions  and  other  errors,  root  mean  square  error  (Erms)  and 
maximum  error  (Emax)  were  calculated.  For  a  mesh  spacing 
of  H  =  0.1,  problem  #1  solutions  proved  to  be  more  accu¬ 
rate.  By  suitably  reducing  the  time  step,  errors  of  less 
than  one  degree  could  be  realized.  No  matter  how  small  a 
time  step  was  used,  solutions  to  problem  #2  were,  at  best, 
accurate  to  within  only  three  degrees.  These  deviations 
occurred  primarily  along  the  grid  boundaries.  An  attempt 
to  reduce  this  error  by  halving  the  mesh  size  resulted  in 
even  larger  error.  This  phenomenon  may  occur  when  large 
numbers  of  arithmetic  operations  cause  the  accumulation  of 
round-off  error  to  become  significant.  Thus,  to  examine 
the  effect  of  mesh  size  on  truncation  error,  additional 
solutions  were  found  for  a  grid  spacing  of  H  =  0.2. 

Root  Mean  Square  Error 

In  both  problems,  the  exponential  method  demon¬ 
strated  superior  accuracy  for  time  steps  greater  than 
0.005  seconds.  This  was  not  surprising  since  this  method 

possesses  the  least  truncation  error  at  small  time  steps 
2  2 

(H  +  At  )  while  enjoying  accuracy  similar  to  the  fully 

2 

implicit  at  large  time  steps  (H  +  At) . 

The  DF  and  SV  schemes  are  comparable  to  the  EXP 
when  time  steps  less  than  0.001  seconds  are  used.  More 
often  than  not,  all  methods  outperform  the  fully  implicit 
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method  for  small  time  steps.  The  DF  scheme  is  superior  to 
the  SV  for  time  steps  below  0.001  seconds.  Above  this 
value,  the  SV  scheme  is  more  accurate.  Dufort  and  Frankel 
point  out  an  optimum  value  of  time  step  for  their  differ¬ 
ence  method.  For  fixed  H,  truncation  error  will  be  mini- 

2  — 

mized  when  a  time  step  of  it  =  H  // 12  is  utilized  (3:142) . 
For  the  mesh  considered  (H  =  0.1),  this  corresponds  to  a 
time  step  of  approximately  0.003  seconds.  Figure  22  sug¬ 
gests  that  minimum  Erms  exist  for  the  DF  method  when  time 
steps  of  0.001  seconds  are  used. 
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Maximum  Error 

Maximum  error,  calculated  at  t  =  0.1  and  t  =  0.3 
seconds,  is  displayed  in  Figures  25-32.  The  effect  of 
doubling  the  mesh  size  can  be  seen  in  Figures  26,  28,  30, 
and  32.  Again,  because  the  exponential  scheme  does  not 
oscillate,  its  performance  is  superior  for  values  of  R 
greater  than  0.5. 

The  DF  scheme  produced  good  results  for  problem  #1 
(Emax  <  one  degree)  when  small  values  of  R  were  used. 
Generally,  the  DF  was  more  accurate  than  the  SV  for  values 
of  R  less  than  0.25.  For  values  of  R  between  0.25  and  0.5, 
the  SV  method  was  more  accurate.  Both  methods  become 
virtually  unusable  for  R  greater  than  1.0  where  solutions 
are  found  to  be  between  10  and  100  degrees  in  error.  Even 
though  the  exponential  solution  retained  superior  accuracy 
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for  large  R,  only  those  solutions  for  R  less  than  1.0  are 
accurate  enough  for  most  engineering  applications. 

An  increase  in  mesh  size  should  result  in  a  larger 
truncation  error  for  all  methods.  For  the  same  R,'  solu¬ 
tions  found  with  H  *  0.1  yield  more  accurate  results  in 
all  cases.  The  effect  of  doubling  the  mesh  size  while  the 

time  step  is  held  constant  can  be  easily  demonstrated.  The 

2  2 

DF  method  possesses  trunction  error  of  order  H  +  At 
2 

+  (At/H)  .  For  t  =  0.0005  and  H  =  0.1  an  error  of 
1.003E-02  is  predicted.  When  H  is  increased  to  0.2,  an 
expected  error  of  4.0007E-02  results.  The  ratio  of  pre¬ 
dicted  errors  is  then  (4 . 0007E-02 ) / (1 . 003E-02 )  or  3.99. 
Figure  27  shows  an  error  of  0.047  for  the  DF  scheme  when 
R  =  0.05  (t  =  0.0005,  H  =  0.1).  The  error  corresponding 
to  this  time  step  on  Figure  28  occurs  at  R  *  0.0125 
(t  =  0.0005,  H  =  0.2)  and  has  a  value  of  0.207.  The  ratio 
of  these  two  errors  is  (0.207/0.047)  or,  approximately  4.40. 

The  effect  of  oscillatory  behavior  on  solution  accu¬ 
racy  is  readily  apparent  from  Figures  29  and  31.  The  SV 
and  DF  solutions  are  reasonably  accurate  (Emax  <  5  degrees) 
until  R  becomes  greater  than  0.5.  Above  R  =  1.0,  oscilla¬ 
tory  behavior  of  both  methods  results  in  significant  error. 

The  performance  of  the  fully  implicit  method  was 
about  what  was  expected.  It  was  superior  to  the  DF  and  SV 
methods  for  large  time  steps.  This  is  because  the  trunca¬ 
tion  error  for  the  DF  and  SV  methods  is  of  order 


H2  +  At2  +  (At/H)2  and  order  H2  +  At  for  the  fully  implicit 
method  (10:189).  This  also  explains  why  the  DF  and  SV 
schemes  are  more  accurate  than  the  fully  implicit  method 
at  smaller  time  steps.  The  exponential  method  was  superior 
to  the  fully  implicit  scheme  for  small  time  steps  and 
roughly  comparable  at  large  time  steps.  Again,  if  the 
truncation  error  of  the  exponential  scheme  is  approxi¬ 
mated  by  that  of  the  Crank-Nicolson  at  small  R,  and  that 
of  the  fully  implicit  at  large  R,  this  behavior  is  expected. 

Time  Required  for  Computations 

All  numerical  solutions  were  calculated  using  a 
24  bit,  Harris  800  computer.  To  compare  the  amount  of 
central  processor  unit  (CPU)  time  required  for  each  method, 
problem  #1  was  solved  using  the  largest  time  step  for  which 
a  maximum  error  of  less  than  one  degree  resulted.  This 
method  of  comparison  is  identical  to  one  used  in  a  previous 
paper  where  the  performance  of  several  finite  difference 
methods  was  also  of  interest  (6:52). 

Of  the  explicit  methods,  the  SV  was  the  fastest. 

The  choice  of  the  fully  implicit  method  to  bootstrap  the 
DF  scheme  to  the  first  time  level  caused  it  to  be  slower 
than  the  SV.  The  EXP  was  slowest  since  a  matrix  equation 
of  order  N,  the  number  of  nodes,  had  to  be  solved  for  each 


iteration . 


fi 


R 

- 


a- 
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For  the  coarser  mesh.  Figure  33  shews  all  methods 
to  be  roughly  the  same.  Both  the  DF  and  SV  required  30 
iterations  to  achieve  the  desired  level  of  accuracy.  The 
EXP  required  only  15  iterations  due  to  its  superior  accu¬ 
racy  at  larger  time  steps. 

When  the  mesh  size  was  halved,  an  obvious  speed 
advantage  resulted  for  the  SV  (Figure  34)  .  The  larger 
matrix  equations  that  resulted  caused  the  EXP  and  DF  methods 
to  be  considerably  slower.  Even  though  the  EXP  again 
required  half  as  many  iterations,  the  SV  scheme  proved 
faster. 

It  is  realized  that  the  grid  sizes  used  here  are 
somewhat  coarse.  But  as  the  mesh  is  refined,  the  SV 
should  increase  its  speed  advantage.  Of  course,  if  R=l/4 
is  used  to  calculate  the  first  time  level  of  the  DF,  the 
method  reduces  to  a  two  level  formula  and  no  additional 
scheme  is  needed  for  bootstrapping.  Interestingly,  the 
apparent  advantage  of  utilizing  larger  time  steps  did  not 
result  in  a  computationally  faster  algorithm  for  the  EXP. 
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Conclusions  and  Recommendations 


Conclusions 

Numerical  oscillations  resulted  in  both  the  DF  and 
SV  schemes  for  the  values  of  R  predicted.  In  both  cases, 
stability  was  maintained  since  the  oscillations  decreased 
with  time.  For  the  same  R,  the  SV  generally  experienced 
less  violent  oscillations.  The  oscillations  were  usually 
large  enough  to  destroy  the  solution  in  both  methods  with 
R  >  1.0.  As  predicted,  the  EXP  did  not  oscillate  for  any 
size  time  step. 

In  problem  #1,  at  early  times  (t  =  0.1  secs),  the 
EXP  retained  good  accuracy  (Emax  <  1  degree)  for  values 
of  R  up  to  0.1.  The  scheme  displayed  marginal  accuracy 
(Emax  <  10  degrees)  for  values  of  R  between  0.1  and  1.0. 

The  exponential  was  more  accurate  than  the  fully  implicit 
method  for  small  time  steps.  For  large  R  (R  >  1.0),  the 
fully  implicit  and  exponential  methods  were  roughly  the 
same.  For  small  R,  the  explicit  methods  and  EXP  were  about 
equal.  The  DF  demonstrated  good  accuracy  for  R  <  0.25. 
Above  this  value,  the  SV  was  consistently  the  more  accu¬ 
rate. 

At  later  times,  the  explicit  solutions  became  more 
accurate  as  the  oscillations  decayed.  The  SV  and  DF  were 
now  accurate  to  within  one  degree  for  values  of  R  up  to 


0.5  (t  =  0.3  secs).  Above  R  =  0.5,  the  SV  was  more  accu¬ 
rate.  In  problem  #1,  the  EXP  retained  good  accuracy  for  R 
up  to  1.0.  Although  the  EXP  was  the  most  accurate  above 
R  =  1.0,  errors  greater  than  10  degrees  still  existed. 

Of  the  explicit  methods,  the  SV  was  the  fastest 
computationally.  The  use  of  the  fully  implicit  method  to 
bootstrap  the  DF  scheme,  caused  it  to  be  fairly  slow  for  an 
explicit  method.  Even  though  the  exponential  method  could 
use  a  larger  time  step  without  sacrificing  accuracy,  the 
burden  of  solving  simultaneous  equations  by  Gaussian 
elimination  caused  it  to  require  the  most  CPU  time. 

Recommendat  cns 

As  a  result  of  this  study,  the  following  are  pro¬ 
posed: 

1.  Generalize  Myer's  matrix  approach  to  oscilla¬ 
tory  behavior  to  include  all  nodes. 

2.  Investigate  the  use  of  the  probabilistic  method 
to  determine  the  stability  condition  of  implicit  finite 
difference  schemes. 

3.  Extend  the  exponential  method  to  three  dimen¬ 
sions  and  check  for  oscillatory  solutions. 

4 .  Apply  the  methods  studied  here  to  problems 
having  derivative  or  mixed  boundary  conditions. 

5 .  Compare  the  accuracy  of  the  exponential  method 
with  that  of  the  variable-weighted  implicit  method  (7:161) 
for  both  small  and  large  time  steps. 
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6 .  Determine  if  the  truncation  error  of  the  exponen¬ 


tial  scheme  can  be  expressed  in  terms  of  R.  If  so,  an 
optimum  value  of  time  step  may  be  realized  resulting  in 
greater  accuracy. 
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Appendix  A:  Derivation  of  the  Duf ort-Frankel 

Finite  Difference  Method 


To  derive  the  Duf ort-Frankel  method  in  two  dimen¬ 


sions,  first  consider  Richardson's  explicit  approximation 


(7:238) : 


T ( i . i . n+1)  -  T(i.i.n-l)  _  T(i+l,j,n)  -  2T(i,j,n)  +  T(i-l,j,nj 
M  AX2 


T ( i, j+1 , n)  ~  2T(i,i,n)  +  T(i,j-l,n) 


This  formulation  is  known  to  be  unconditionally 


unstable  (7:238).  If  the  T(i,j,n)  term  is  replaced  by  (1/2) 
(T ( i, j ,n+l)  +  T (i, j ,n-l) )  the  unconditionally  stable  Duf ort- 


Frankel  scheme  results: 


T(i,j,n+1) (1+4R)  =  2R[T(i+l, j,n)  +  T(i-l,j,n) 


+  T(i, j+l,n)  +  T(i, j-l,n) ] 


+  ( 1-4R)  T ( i, j , n-1) 


where 


R  =  At/H‘ 


Ax  =  Ay  =  H 


R 

& 
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Appendix  B :  Matrix  Analysis  of  Oscillations  in  the 

Exponential  Difference  Me thod 

Myer's  matrix  approach  to  predicting  oscillatory 
behavior  will  be  applied  to  the  exponential  scheme  in  two 
dimensions.  In  addition,  a  four  nodal  point  grid  and 
homogeneous  Dirichlet  boundary  conditions  are  assumed. 

To  begin,  a  difference  equation  is  written  for  each 
node  of  the  grid.  The  resulting  system  of  four  equations 
can  then  be  written  as: 


A  u(n+1)  =  B  H(n) 


K 


s 


M  * 

V\ 

*  ft 


Here,  A  and  B  represent  coefficient  matrices,  while  u^n+1^ 
and  u(n)  are  column  matrices  of  temperatures  at  times  (n+1) 
and  (n)  respectively.  Further,  the  eigenvalues  of  A  ^  B 
will  determine  the  oscillatory  behavior  of  the  system  (8:285) 
and  are  found  from  the  relation  (8:285) 

det (A-1  B  -  XI)  =  0 

Myers  shows  that  this  is  equivalent  to  the  expres¬ 


sion  (8:285) 

det  (B  -  XA)  =  0 

Hence,  the  eigenvalues  of  (B  -  XA)  determine  the 
oscillatory  behavior  of  the  system.  If  Equation  (3)  is 
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applied  to  a  four  nodal  point  grid  with  homogeneous 
Dirichlet  boundary  conditions,  the  following  matrix  results 
for  (B  -  XA) : 


i 

1 


;v 


2 

$ 

i 

8 

R 


X  -  \Y 

C  +  XD 

C  +  XD 

0 

C  +  XD 

X  -  XY 

0 

C  +  XD 

C  +  XD 

0 

X  -  XY 

C  +  XD 

0 

C  +  XD 

C  +  XD 

X  -  XY 

where 

X  =  1  -  4R ( 1-f ) 

Y  =  1  +  4Rf 
C  *  R ( 1-f ) 

D  =  Rf 

The  determinant  can  be  simplified  using  the  follow¬ 
ing  techniques  (4) : 

Subtracting  the  last  row  from  the  first  yields 


X  -  XY 

0 

0 

-(X  -  XY) 

C  +  XD 

X  -  XY 

0 

C  +  XD 

C  +  XD 

0 

X  -  XY 

C  +  XD 

0 

C  +  XD 

C  +  XD 

X  -  XY 
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Subtracting  the  third  row  from  the  second  and  factoring 
results  in 


(X  -  AY) 


C  +  AD 


X  -  AY 


C  +  AD 


C  +  AD 


C  +  AD 


X  -  AY 


Adding  the  first  column  to  the  last  yields 


(X  -  AY) 


C  +  AD 


X  -  AY 


2 (C  +  AD) 


C  +  AD 


C  +  AD 


X  -  AY 


Developing  the  determinant  about  the  first  row 


results  in 


(X  -  AY) 


X  -  AY 


C  +  AD 


2 (C  +  AD) 


c  +  ad  X  -  ay 


Adding  the  first  column  to  the  second  and  again  developing 
about  the  first  row  yields 


X  -  XY  2  (C  +  XD) 

(X  -  XY)2 

2  (C  +  XD)  X  -  XY 

For  the  value  of  the  determinant  to  be  zero, 
the  following  relations  must  be  true 

(X  -  XY)  2  =  0 

(X  -  XY)  ■  +  2(C  +  XD) 

(X  -  XY)  =  -  2 (C  +  XD) 

The  first  relation  produces  two  identical  eigenvalues, 

(X/Y)  or: 

(l-4R(l-f ) ) / (l+4Rf)  =  exp  (-4R) 

The  second  and  third  relations  respectively  yield 

X3  -  1  -  (6R/ (l+6Rf) ) 

and  X4  =  1  -  (2R/ (l-2Rf ) ) 

It  is  easily  verified  that  all  eigenvalues  are  less 
than  one  in  modulus  and  that  X4  is  the  spectral  radius. 

For  all  values  of  R,  X4  is  positive  and  less  than  unity, 
hence,  no  oscillations  are  predicted. 
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